/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
//
//  Xianguo Lu <lu@physi.uni-heidelberg.de>
//

#include <TFile.h>
#include <TH2D.h>
#include <TVector3.h>

#include "AliTPCseed.h"
#include "AliTRDtrackV1.h"

#include "AliESDfriend.h"

#include "AliCosmicTracker.h"
#include "AliESDCosmicTrack.h"
#include "AliTPCCosmicUtils.h"
#include "AliTRDdEdxBaseUtils.h"

#include "AliAnalysisTaskdEdxCosmic.h"

AliAnalysisTaskdEdxCosmic::AliAnalysisTaskdEdxCosmic(const TString chunk):
  AliAnalysisTaskSE("dEdxCosmic")
  , fStream(0x0)

  , fESDEvent(0x0)
  , fESDfriend(0x0)

  , fChunk(chunk)

  , fhtrig(0x0)
  , fhntrdv1(0x0)
  , fhnpair(0x0)
  , fhreuse(0x0)

  , fhminphi(0x0)
  , fhnfillevent(0x0)
  , fhnfilltrk(0x0)

  , fHist(0x0)
{
  //
  //
  //
  fhtrig = new TH1I("a0_htrig", "", 40,0,40);
  fhntrdv1 = new TH1I("a1_hntrdv1","",20,0,20);
  fhnpair = new TH2I("a25_hnpair","",50, 0, 50, 100,0,100);
  fhreuse = new TH1I("a26_hreuse","",2,0,2);

  fhminphi = new TH1D("a5_minphi","",90,0,90);
  fhnfillevent = new TH1I("a6_nfillevent","",20, 0, 20);
  fhnfilltrk = new TH1I("a7_nfilltrk","",20, -10, 10);

  fHist = new TList;
  fHist->SetOwner();//------>>> Very important!!!!!!!
  
  fHist->Add(fhtrig);
  fHist->Add(fhntrdv1);
  fHist->Add(fhnpair);
  fHist->Add(fhreuse);
  fHist->Add(fhminphi);
  fHist->Add(fhnfillevent);
  fHist->Add(fhnfilltrk);

  fStream = new TTreeSRedirector(Form("TaskdEdxCosmicTree_%s.root", fChunk.Data()));
}

//____________________________________________________________
AliAnalysisTaskdEdxCosmic::~AliAnalysisTaskdEdxCosmic(){
  //
  // Destructor
  //
  delete fStream;
  delete fESDEvent;
  delete fESDfriend;

  fHist->Clear();
  delete fHist;
}

void AliAnalysisTaskdEdxCosmic::UserExec(Option_t *)
{
  //
  //Execute
  // 
  static Int_t ievt = 0;
  //if(ievt%100==0) printf("iEvent: %d\n", ievt);
  ievt++;

  if(!AliTPCCosmicUtils::GetESD(fESDEvent, fESDfriend)){
    printf("xlulog no esd!!\n"); exit(1);
  }

  Int_t mag = -999, btrig = -999, run=-999, spe=-999;
  if(!CutESD(mag, btrig, run, spe))
    return;

  //const TClonesArray *stackCosmic = AliCosmicTracker::FindCosmic(fESDEvent, kFALSE);

  //Int_t npair = stackCosmic->GetEntries();
  Int_t npair = fESDEvent->GetNumberOfCosmicTracks();

  fhnpair->Fill(fESDEvent->GetNumberOfTracks(), npair);


  for(Int_t ip=0; ip<npair; ip++){
    //const AliESDCosmicTrack * esdcos = (AliESDCosmicTrack *) stackCosmic->At(ip);
    const AliESDCosmicTrack * esdcos = fESDEvent->GetCosmicTrack(ip);

    fhreuse->Fill(esdcos->IsReuse());
    if(esdcos->IsReuse())
      continue;

    AliESDtrack *trks[2]={fESDEvent->GetTrack(esdcos->GetESDUpperTrackIndex()), 
                                fESDEvent->GetTrack(esdcos->GetESDLowerTrackIndex())};

    const AliExternalTrackParam * params[]={esdcos, esdcos->GetLowerPartner()};

    if(!CutMinPhi(params))
      continue;

    //----------------------- filling --------------------------
    Int_t nfillevent=0;
    for(Int_t itrk=0; itrk<2; itrk++){
      AliESDtrack * esdtrack = trks[itrk];

      AliTPCseed * tpcseed = AliTPCCosmicUtils::GetTPCseed(esdtrack);
      if(!tpcseed){
        fhnfilltrk->Fill(0);
        continue;
      }

      AliTRDtrackV1 * trdv1 = AliTRDdEdxBaseUtils::GetTRDtrackV1(esdtrack);
      if(!trdv1){
        fhnfilltrk->Fill(-1);
        continue;
      }

      const Int_t ntrklet = trdv1->GetNumberOfTracklets();
      if(!ntrklet){
        fhnfilltrk->Fill(-2);
        continue;
      }

      const Bool_t ksame = AliTRDdEdxBaseUtils::IsInSameStack(trdv1);
      if(!ksame){
        fhnfilltrk->Fill(-3);
        continue;
      }

      Int_t istack = -999, isector = -999;
      Double_t mom = -999;
      AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(trdv1, isector, istack, mom);
      if(istack<0){
        printf("xlulog error AliTRDdEdxUtils::GetFirstSectorStackMomentum() fail!! %d %d %f\n", isector, istack, mom); exit(1);
      }

      Double_t ltbg = params[itrk]->P();
      //anchor at bg=25
      if(ltbg<2){
        ltbg = mom;
      }
      ltbg = TMath::Log10(ltbg/AliTPCCosmicUtils::Mass());

      Double_t tpcsig  = esdtrack->GetTPCsignal();

      Int_t pid = AliPID::kMuon;

      Int_t charge = params[itrk]->Charge();

      (*fStream)<<"tree"<<
        "npair=" << npair<<

        "run="<<run<<
        "itrk="<<itrk<<
        "trig="<<btrig<<
        "mag="<<mag<<
        "charge="<<charge<<
        "spe="<<spe<<
        "pid="<<pid<<

        "isector="<<isector<<
        "istack="<<istack<<

        "ltbg="<<ltbg<<

        "esdtrack="<<esdtrack<<
        "trdv1="<<trdv1<<
        "tpcseed="<<tpcseed<<
        "tpcsig="<<tpcsig<<

        "\n";

      nfillevent++;
      fhnfilltrk->Fill(1); 
    } 
    fhnfillevent->Fill(nfillevent);
  }

  //delete stackCosmic;
}

void AliAnalysisTaskdEdxCosmic::Terminate(Option_t *) 
{
  //
  //do when Terminate
  //
  TFile *fout = fStream->GetFile();
  if(fout){
    fout->cd();
    fout->Write();
    fout->Close();
  }

  TFile *fh = new TFile(Form("TaskdEdxCosmicHist_%s.root", fChunk.Data()), "recreate");
  fHist->Write();
  fh->Save();
  fh->Close();
  delete fh;
}

//========================================================================================
//========================================================================================

Bool_t AliAnalysisTaskdEdxCosmic::CutESD(Int_t & mag, Int_t &btrig, Int_t &run, Int_t &spe)
{
  //
  //
  //

  //AliESDHeader::fEventSpecie Reconstruction event specie (1-default,2-lowM,4-highM,8-cosmic,16-cal)
  spe = fESDEvent->GetEventSpecie();

  run = fESDEvent->GetRunNumber();
  mag = AliTPCCosmicUtils::GetBField(fESDEvent);

  btrig = AliTPCCosmicUtils::GetTrigger(fESDEvent)&(AliTPCCosmicUtils::k0HWU+AliTPCCosmicUtils::kTRDCO2);
  fhtrig->Fill(btrig);
  if(!btrig){
    return kFALSE;
  }

  const Int_t ntrdv1= AliTRDdEdxBaseUtils::GetNtracklet(fESDEvent);
  fhntrdv1->Fill(ntrdv1);
  if(!ntrdv1){
    return kFALSE;
  }

  return kTRUE;
}

Bool_t AliAnalysisTaskdEdxCosmic::CutMinPhi(const AliExternalTrackParam *params[])
{
  const Double_t minphi = AliTPCCosmicUtils::GetMinPhi(params);
  const Double_t mpcut = 10;
  fhminphi->Fill(minphi);
  if(minphi<mpcut)
    return kFALSE;
  else
    return kTRUE;
}
